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We introduce a 'virtual-move' Monte Carlo (VMMC) algorithm for systems of pairwise-interacting 
particles. This algorithm facilitates the simulation of particles possessing attractions of short range 
and arbitrary strength and geometry, an important realization being self-assembling particles en- 
dowed with strong, short-ranged and angularly specific ('patchy') attractions. Standard Monte 
Carlo techniques employ sequential updates of particles and can suffer from low acceptance rates 
when attractions are strong. In this event, collective motion can be strongly suppressed. Our al- 
gorithm avoids this problem by proposing simultaneous moves of collections (clusters) of particles 
according to gradients of interaction energies. One particle first executes a 'virtual' trial move. We 
determine which of its neighbours move in a similar fashion by calculating individual bond ener- 
gies before and after the proposed move. We iterate this procedure and update simultaneously the 
positions of all affected particles. Particles move according to an approximation of realistic dy- 
namics without requiring the explicit computation of forces, and without the step size restrictions 
required when integrating equations of motion. We employ a size- and shape-dependent damping of 
cluster movements, motivated by collective hydrodynamic effects neglected in simple implementa- 
tions of Brownian dynamics. We discuss the virtual-move algorithm in the context of other Monte 
Carlo cluster-move schemes, and demonstrate its utility by applying it to a model of biological 
self-assembly. 



I. INTRODUCTION 

A standard Monte Carlo simulation of a system of par- 
ticles consists of a sequence of moves of individual parti- 
cles. If these moves are proposed and accepted according 
to the principle of detailed balance (or balance [H or 
superdetailed balance [3j), the system will eventually re- 
lax to thermal equilibrium 0]. For some models, this 
approach even provides an approximation to the dynam- 
ics that the corresponding physical system would exe- 
cute [SJ|5J|7]. However, such approximations break down 
when relaxation requires the movement of particles in 
concert, particularly in the presence of very strong inter- 
actions. This is the case, for example, when small ions 
of opposite charge bind together tightly [Sj. 

In such situations, the suppression of collective modes 
of motion exhibited by standard Monte Carlo algorithms 
affects both intra-cluster relaxation and whole-cluster 
diffusion. As an example of the latter we show in Figure [l] 
four particles equipped with strong pairwise interactions. 
Only two particles (j and k) are close enough to interact. 
An 'overdamped' (Langevin) dynamics would see particle 
positions evolve according to both deterministic forces, 
derived from pairwise interactions, and random buffeting 
forces designed to model solvent fluctuations. The ran- 
dom forces induce collective diffusion of isolated clusters 
(such as the dimer jk). Molecular dynamics simulations 
naturally accommodate collective motion by updating si- 
multaneously the position of every particle, according to 
both deterministic and random forces. A standard Monte 
Carlo simulation, however, employs sequential updates of 
individual particle positions. The effects of the determin- 



istic and random buffeting forces are then inextricably 
linked. Potential energy gradients dictate that moves of 
one particle relative to another (e.g. the displacement 
that changes configuration a to configuration b) are sup- 
pressed by the exponential of the change in interaction 
energy. If this change is but a few multiples of k^T, 
such a move is unlikely to be accepted. Since the ran- 
dom buffeting force is modeled by the random motion 
of individual particles, rejecting such motion has the un- 
desirable effect of annulling the instantaneous buffeting 
force experienced by that particle. This kinetic 'trap' 
leads to a suppression of collective diffusion, and hence 
to an unphysical dynamics. 

In a similar manner, internal cluster relaxation in- 
volving collective modes of motion is under-represented 
by making sequential moves of single particles. A re- 
alistic description of collective structural rearrangement 
is necessary to model many phase transitions, such as 
the microphase separation of colloids at low tempera- 
ture, and examples of self-assembly, e.g. of large pro- 
teins called chaperonins [THl IZJ [TH]; various nanoparti- 
cles [inilSnjIlI]; and virus capsids P^E5]. 

For sufficiently small displacements, single-particle 
moves can in principle retain acceptance rates large 
enough to permit some collective modes of relaxation. 
However, when interactions are not only strong but also 
short-ranged, the displacement scale required to ensure 
acceptable rates of collective motion is invariably so small 
that configuration space cannot be explored in a reason- 
able time. 

Compounding these difficulties in modeling a natural 
dynamics is the fact that even in the limit of unit ac- 
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clusters with greater flexibility can lead to efficent collec- 
tive re-arrangement. However, while equilibration can in 
this manner be facilitated, identifying and moving clus- 
ters according to properties of the current configuration 
does not restore a physical dynamics: particles do not 
explore local gradients of potential energy in a realistic 
manner. We discuss this situation in more detail in Ap- 
pendices A and B. 



FIG. 1: An illustration of the difficulty encountered by stan- 
dard Monte Carlo methods in the presence of very strong, 
short-ranged interactions. Here only particles j and k inter- 
act, and do so strongly. Single-particle Monte Carlo schemes 
generate a diffusion of the dimer jk via relative moves of j and 
k. Because such moves (e.g. a ^ b) are suppressed by the 
exponential of the resulting energy change, collective modes 
of motion are under-represented for any acceptance rate less 
than unity. Moving particles j and k simultaneously (a c) 
can restore collective motion. We shall discuss how this can 
be done in order to approximate a realistic dynamics. 



ceptance rate the diffusion associated with single-particle 
Monte Carlo protocols represents only a rough approx- 
imation of physical motion. Single-particle moves (and 
indeed standard implementations of Brownian dynam- 
ics) induce motion of an isolated cluster with a trans- 
lational diffusion constant that scales as the reciprocal 
of the number of constituents of the cluster, and a rota- 
tional diffusion constant proportional to the reciprocal of 
the cluster's moment of inertia [9 . By contrast. Stokes' 
law for a body moving through a viscous medium implies 
translational and rotational diffusion properties depend- 
ing on the first and third powers, respectively, of the 
cluster's characteristic radius. 

These limitations can be overcome in principle by 
proposing simultaneous moves of clusters of particles (e.g. 
the move that takes state a to state c in Figure Ex- 
isting cluster algorithms [lOj [TTJ [121 HSl HI] have been 
used with great success to simulate many complex sys- 
tems [TS]. Such algorithms define a cluster to be moved 
by recursively linking particle pairs with a given prob- 
ability. For those algorithms that effect local moves of 
particles, the simultaneous displacement of collections of 
particles allows the restoration of diffusive modes of mo- 
tion suppressed by single-particle schemes. 

However, such algorithms restore diffusive motion by 
identifying and moving clusters based on the properties 
of particles (energy or proximity) in their current config- 
uration, without regard for the changes induced by the 
proposed move. As a consequence, these algorithms tend 
to under-represent internal collective motion: particles 
interacting strongly (or particles in close proximity) will 
be displaced collectively, and will therefore not move rel- 
ative to each other. This leads to the development of 
severe kinetic traps. In Appendix A we discuss how gen- 
eralizing these algorithms to permit the identification of 



To permit both realistic diffusive motion and collective 
internal re-arrangement we propose defining and mov- 
ing clusters on the basis of gradients in potential energy. 
In this paper we introduce a 'virtual-move' Monte Carlo 
(VMMC) scheme designed to approximate a realistic dy- 
namics for particles possessing strong, short-ranged in- 
teractions. In our scheme, one particle first executes a 
'virtual' trial move. We determine which of its neigh- 
bours move in a similar fashion by calculating bond en- 
ergies before and after the proposed move. We iterate 
this procedure from all 'recruited' neighbours and stop 
when no more particles show a tendency to move. We 
then update simultaneously the positions of all affected 
particles. 



The VMMC scheme is designed to facilitate the sim- 
ulation of components that spontaneously self-assemble. 
Such components typically possess pairwise attractions 
of arbitrary strength, possibly very short range, and a 
high degree of angular specificity ('patchiness'). The al- 
gorithm executes collective updates of particles in a man- 
ner designed to approximate realistic particle motion. It 
permits a basic particle displacement scale larger than 
can be used, for example, when integrating equations of 
motion. Furthermore, cluster moves may be performed 
so as to respect Stokes' law, offering the possibility of ex- 
ceeding the dynamical realism of simple implementations 
of Brownian dynamics. 



This paper is organised as follows. In Section |TT] 
we show that by proposing collective moves of parti- 
cles based upon individual bond energy gradients one 
can evolve in an approximately realistic (and compu- 
tationally efficient) way a system of strongly pairwise- 
interacting particles. In Section |III| we generalize the 
scheme to permit the use of distinct real and virtual 
moves, thereby allowing precise control of relative ro- 
tational and translational motion. We show that in 
this way we can mimic the types of relaxation observed 
in Brownian dynamics simulations of isotropic Lennard- 

we ap- 



Jones particles in two dimensions. In Section IV 



ply the VMMC algorithm to a system of self-assembling 
components designed to model the aggregation of protein 
complexes called chaperonins. We conclude in Section [V] 
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II. A 'VIRTUAL-MOVE' MONTE CARLO 
CLUSTER ALGORITHM 

A. Making collective moves 

In this section we introduce our central result, a 
virtual-move Monte Carlo algorithm. This algorithm is 
designed to approximate physical dynamics by moving 
interacting bodies either individually or in concert ac- 
cording to individual bond energy gradients. This proce- 
dure may be regarded as a particle-based adaptation of 
the Wolff cluster algorithm [T3] , and may be used on- or 
off-lattice. It is similar in some respects to the algorithm 
described in Ref. [lOj, but differs in that here proposed 
moves are chosen to correspond as closely as possible to 
realistic particle displacements. 

As discussed in the introduction, making sequential 
moves of single particles leads to a suppression of col- 
lective modes of motion. To correct this deficiency it is 
necessary to effect moves of particles in concert. We shall 
explain in this subsection how this can be done in gen- 
eral. In the following subsection we describe the specifics 
of the virtual- move scheme. In brief, we effect cluster dis- 
placements by first proposing a move of a single particle. 
If the change in energy of interaction between that parti- 
cle and its neighbours is unfavourable, those neighbours 
are 'recruited', adopting the move of the first particle as 
their own. We continue this recruitment until no further 
particles show a tendency to move. Before describing 
this scheme in detail we shall set up our notation and 
nomenclature. 

We consider in d dimensions a system of N particles 
having radius Rq and interacting by way of a short- 
ranged, pairwise potential e^. In Figure [2] we show a 
typical cluster move (and its reverse) from state fi to 
state ly. We define two particles to be contiguous if their 
centres are separated by a distance less than or equal 
to Tc, the interaction range of the potential. We define 
a physical cluster (abbreviated simply to 'cluster') as a 
group of contiguous particles. We shall use an iterative 
linking procedure, described below, to select a group of 
particles, C. This group can range in size from a sin- 
gle particle to any one cluster in the system. We shall 
call this group of particles a pseudocluster. In Figure [2] 
the pseudocluster C is shaded. We perform on C either a 
rigid-body rotation or a translation in an arbitrary direc- 
tion (or, if desired, a linear combination of a translation 
and a rotation). In the figure, one such translation hap- 
pens to bring the pseudocluster into contact with the 
cluster C2. 

To define a pseudocluster C we select from the sys- 
tem a 'seed' particle i. The seed is the first member of 
the pseudocluster. We shall attempt to 'recruit' particles 
as members of the pseudocluster by forming 'links' from 
members of the pseudocluster to neighbouring particles. 
We start by forming a link with probability Pij{fJ- — > i') 
between i and any particle j with which i interacts. The 
form of pij shall be specified later. If this link forms we 
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FIG. 2: An illustrative collective move. We define a pseudo- 
cluster C (shaded particles) using an iterative linking scheme 
(see text). A particular realization of links is denoted by bold 
black lines. Failed links internal to the pseudocluster are not 
shown. We subject the pseudocluster to a rotation or a trans- 
lation (or both). Here a proposed translation moves C from 
contact with cluster Ci in state /j, to contact with cluster C2 in 
state u. The interface between the pseudocluster and its envi- 
ronment in state a is labelled 1^, and is defined by the set of 
all pairwise interactions between white and shaded particles. 

regard j as a member of the pseudocluster. We continue 
the linking procedure until links have been proposed ex- 
actly once from each member m of the pseudocluster to 
all particles with which m interacts. 

Linking particles in this manner results in the pseudo- 
cluster C being selected and moved with probability 

Wgon(M ^ l^) = ^'sccd(M) ^displace (C; ^J, ^ v) ■ (1) 

Here Psoed(M) accounts for the likelihood of choosing 
a seed particle in the pseudocluster C. The factor 
^dispiace(C; ^ ^ v) retums the probability (given a seed 
particle) of building the pseudocluster in state /i and 
moving it such that the proposed new state is v. It de- 
pends on two factors: 

c TZ 

^dispiaco(c; A* ^ ^) = X! n ^ n ~^ 

The first factor in Equation (p| concerns the link-failure 
probabilities qij{fJ. —^i') = l ^ pij{ii v). The product 
O/i^i/ runs over all links which must not form in order to 
move from state ^ to state u; unformed links external to 
the pseudocluster define the interface labeled in Fig- 
ure [2] The second factor in ^ is the probability of form- 
ing links (labelled €) between the constituent monomers 
of C. In general, there are many distinct realizations TZ of 
links (and failed links) leading to the same chosen pseu- 
docluster (denoted by X^k)' ^^'^ product runs over 
links comprising one such realization. Again in general, 
the probability of forming links between particles such 
that the shaded cluster is the chosen pseudocluster de- 
pends on whether one is executing the forward or reverse 
move. 

The concept of detailed balance ensures that a system 
evolves towards equilibrium by requiring that the rates 
W for passing between any two states and v satisfy 

p{^l)w{^l^v) = p{v)w{v^ (a) 
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Here p{a) = Z~^e~^^° is the Boltzmann weight with 
respect to the system's Hamiltonian in state a G {/i, z^}, 
Ha', P = 1/T is the reciprocal temperature (we adopt 
units such that fee = 1); and Z is the partition function. 
Calculating the rates W would require the enumeration 
of all possible ways of linking together monomers in order 
to generate a given pseudocluster. For large systems this 
is not feasible. Instead, we use an iterative procedure 
to identify one particular realization of links and failed 
links, and compute the likelihood that each of these links 
(or failed links) will form when making the reverse move. 
Wc then impose the requirement 



p{fi)W{fi ^ iy\TZ) = p{u)W{v ^ p\n), 



(4) 



where W{p — > v\'R-) is the rate for passing from state p 
to state V, given a realization TZ of links and failed links. 
This is the condition of superdetailed balance [3]. The 
realization TZ must include the direction in which links 
are formed (e.g. from particle i to j, or vice versa). 

We note finally that the rate for passing between states 
is the product of the rates of proposing (generating) and 
of accepting such a move: W = Wgcn x W^acc- 

Combining the results of this section fixes the ratio 
of acceptance rates for forward and reverse moves. We 
choose as the acceptance rate for the move p v 

Wacc {P ^ 

= Q (Tie - nc)V(C) mm i I, — -— c '^^ " 



n^.^^ q,j{v^ P) -pr p„{v^ p) 



(5) 



The factors outside the 'minimum' function are not re- 
quired to enforce balance (and therefore to ensure sam- 
pling the correct equilibrium distribution), but are cho- 
sen to preserve as closely as possible a realistic dynamics. 
The first factor in Equation ([s]) is used in conjunction 
with an 'early rejection' termination of the link-forming 
procedure if we recruit more particles than is 'physical'. 
The function Q returns zero if the number of particles in 
the pseudocluster, nc is greater than a specified cutoff ric, 
and unity otherwise. In order to propose moves of clus- 
ters of arbitrary size with correct frequency it is necessary 
to suppress moves of clusters of size nc > 1 by drawing 
the cutoff Tic from a particular distribution. Otherwise, 
large clusters move in general more frequently than small 
clusters. It is crucial to enforce this frequency-correction 
during the link-forming procedure; we would waste much 
time by simply rejecting moves of large clusters with high 
probability. 

The factor I?(C) < 1 encodes the diffusion properties 
of the pseudocluster C. One of the advantages of the 
Monte Carlo cluster-move framework is that we know in 
advance the size and shape of the aggregate C whose dis- 
placement we are considering. We can then enforce a 
hydrodynamic damping of translations and rotations by 



rejecting moves of clusters with a factor that depends 
on the size and shape of the aggregate. In Section III 
we discuss how to scale collective translations and rota- 
tions instead of simply rejecting collective moves with a 
specified probability. 

The first factor (of the second argument) within the 
braces accounts for the likelihood of picking a given par- 
ticle i as a seed in states p and v. We shall choose the 
seed particle uniformly from any in the system, giving 

Psccd{p) = -Pseed(l^)- 

To derive the final factor of the first line within the 
braces we have used the result that the internal pseu- 
docluster energies in states p and v are identical, and 
so the ratio of Boltzmann weights p{i') / p{p) reduces to 
Q~P{E„-E^) ^ Here is the interfacial energy between 
the cluster and its environment in state a, namely, the 
sum in that state of all pairwise energies between white 
and shaded particles. 

The final line of Equation (|5| contains the link-failure 
{qij) and link- forming {pij) factors. These depend on 
the specific choice of the linking probability pij {p v). 
In the following subsection we discuss in detail one such 
choice. Note that the link- failure products run over failed 
links both internal to and external to the pseudocluster. 
The link-making factor is evaluated for a given realization 
TZ of links. 



B. Making collective moves according to potential 
energy gradients 

With the Monte Carlo cluster-move framework de- 
scribed, we turn to a specific choice for the hnking proba- 
bility Pij (/i — *■ v). To enforce an approximate dynamical 
realism this probability must depend on the proposed 
move connecting the initial state p to the final state v. 
In Appendices A and B we discuss how linking particles 
according only to properties of the initial state leads to 
a dynamics that is not physical. 

Our virtual-move scheme is as follows. We start in 
state p. We select uniformly a pseudocluster 'seed' par- 
ticle i having coordinates (position and orientation) Xi, 
and assign to this seed a move map. This map defines 
a random translation or rotation about an axis through 
the centre of the seed (or, if desired, a linear combination 
of both) . We denote by x'^ the coordinates of i following 
application of this map. We execute a 'virtual' move of 
the seed under the map, and look to those particles {j} 
with which the seed interacts in state p. We link a given 
particle j to the seed with a probability Pij{p v) that 
depends on the energy of the relevant bond before and 
after the move of the seed: 

Ptj{p ^ = ("•£ - nc) 

X 2i;Vax (0, 1 - e''^^(''^")-'5^'(''J")) . (6) 

The energy 

Ec{i, j) = e^j {x[, x'j) = ey (x;, Xj) (7) 



is the energy of the bond ij following a collective virtual 
move of i and j (where each move according to the same 
map) . Because the map defines either a rigid-body rota- 
tion or translation, the particles do not move relative to 
each other, and so this bond energy is identical to that in 
the starting state /i. This gives rise to the second equality 
in Equation (|7]). The term 



(Xj, Xj) 



(8) 



is the bond energy following an individual virtual move 
of i, and no move of j (the bond energy in the proposed 
state i^). 

We define the 'interaction' term j'j^^^ in Equation mi 
to be unity if in state particles i and j are deemed 
to be interacting (and therefore eligible for the linking 
procedure), and zero otherwise. We define interacting 
particle pairs as those whose centres are separated by less 
than the range of the interaction. We shall see later that 
we must account separately, via the overall acceptance 
rate, for certain particle pairs that start or end a move 
in a 'noninteracting' configuration. 

The first factor in Equation ^ enforces an 'early re- 
jection' of the link-forming procedure. Because collective 
moves can in principle be initiated from any particle, 
those particles residing in large clusters have a greater 
chance of changing position than do isolated particles. 
We account for this by suppressing the rate for moves of 
clusters of size nc by a factor of 1 /nc . It is not efficient 
to engineer this suppression via the acceptance rate. A 
priori we do not know how many particles will be as- 
signed to a pseudocluster. We would waste much effort 
by building large pseudoclusters only to reject their move 
with high probability. Instead, we suppress the genera- 
tion rate for pseudocluster construction. For each move 
we draw a cutoff ric from the distribution Q{nc) oc n~^, 
and abort the link formation procedure if the pseudoclus- 
ter size exceeds ric- We then cancel the move, as indicated 
by the prefactor in Equation ([s]). In this way we ensure 
(with reasonable computational efficiency) that particles 
experience positional updates with approximately equal 
frequency, a requirement necessary for dynamical fidelity. 
The early-rejection procedure can be further modified to 
account for the diffusion properties of clusters (see fol- 
lowing section). 

If a link is formed, the linkee particle j is 'recruited' 
to the pseudocluster, and adopts the move map of the 
seed. The linker particle i is then returned to its original 
position. We next perform a 'reverse' virtual move of the 
linker particle and record the probability Pij{iy — > /i) that 
the link ij would form were the linker particle translat- 
ing in the opposite direction (and/or rotating with the 
opposite sense). This factor will be used to enforce bal- 
ance in the manner shown in Equation ([s]); it suppresses, 
for example, collective moves of hard particles with no 
attractive interactions j26] . 

We continue this scheme hierarchically, attempting to 
link (exactly once) each member m of the pseudocluster 
to any unlinked particle with which m interacts. We stop 
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FIG. 3: An illustration of virtual-move Monte Carlo. From 
starting state fi the seed particle i is assigned a virtual move 
map, denoted by a thin arrow (a). We propose a link between 
i and j by calculating the energy difference of the bond ij 
before and after the virtual move of j (b, c). In this example 
a link forms, and j adopts the virtual move map of i; the 
latter is returned to its original position (d). We iterate this 
procedure until no more links form (e), and displace all linked 
particles simultaneously (bold arrows denote real, not virtual, 
moves): (f) — > v. The new configuration u is proposed as the 
final state, and the Monte Carlo acceptance probability is 
evaluated. 



when no further links form. We then update simulta- 
neously the positions of all linked particles, moving the 
pseudocluster as a rigid body, and evaluate the Monte 
Carlo acceptance probability. 

In intuitive terms, with probability Pij{fJ. v) we pro- 
pose a move of i and j in concert, according to a move 
map specifying a change in position of i from a state /i 
to a state v. With the complementary probability we in- 
stead propose a move of i relative to j. This is the key dif- 
ference from a standard single-particle MC scheme: here, 
if the move of particle i relative to j is rejected, i and j 
are moved in concert. 

The scheme is best explained using an example, which 
we show in Figure [3j We consider five particles, i, j, 
fc, / and m, endowed with short-ranged, orientation- 
dependent pairwise interactions (particle orientations not 
shown) . Bonds ij, ik and jk are strong, possessing a large 
negative energy e, while bond kl is weak, possessing an 
energy of interaction (5 « 0. There is no interaction be- 
tween m and any other particle. 

We first choose a seed particle, say i. To i we assign a 
move map, denoted by a thin arrow (a). The map here 
defines a rightward translation. We shall denote the state 
of particle a by — {^a, Sq,} (which enodes the position 
Va and orientation Sq, of the particle). After executing 
its virtual move, the state of particle a is x^. 

To determine whether i moves individually or in con- 
cert with other particles, we form links between i and 
each particle with a probability given by Equation Q . In 



our example the linking procedure unfolds as follows. We 
propose a link between i and j , and so consider only these 
particles. We calculate the initial energy of the bond ij, 
Ec, and find it to be large and negative Ec = e (b). 
We move i under its map and calculate the new energy 
of the bond (c). Because attractions are short-ranged 
we find this energy to be zero. We link i and j with 
probability — > ^) = 1 — e^*^ w 1. In our example 
this link is accepted, and so j adopts the move map of 
i. The latter is returned to its original position (d). We 
then execute a 'reverse' virtual move of i (not shown), 
and record the reverse linking probability Pij{i^ — > fJ-). 
We next propose a link between j and k, and find that 
it is also accepted (e). Again, we record the reverse link- 
ing probability. We continue the link-forming procedure 
by testing the interaction kl; this weak bond is linked 
with probability pki ~ S/k^T w 0. In our example this 
link is not formed. The virtual-move linking procedure 
is now finished. The final virtual moves are adopted as 
real moves (thin arrows become bold arrows, panel (f)) 
and all particles in the chosen pseudocluster {i,j,k) are 
displaced simultaneously, according to their map (defin- 
ing a new state, v). We propose this state v as the final 
state, and evaluate the Monte Carlo acceptance factor. 

The 'reverse' virtual move described above is used to 
ensure superdetailed balance by suppressing the likeli- 
hood of a move /i — > i/ by the product over all links of 
terms of the form Piji^v fi)/pij{fi — s- v), as shown in 
Equation ^ . The necessity of this procedure can be un- 
derstood in intuitive terms by referring to Figure |4j In 
this example, the forward move f is initiated by displac- 
ing the seed i rightwards. This move causes i to overlap 
particles j and k. If these overlaps are 'hard' (infinitely 
unfavourable energetically), links ij and ik form with cer- 
tainty; we do not need to propose a link between k and 
j. The chosen pseudocluster C is the trimer ijk, and the 
proposed final state follows by displacing C rightwards. 

But now we encounter a problem generating the re- 
verse move. We cannot do so by displacing, for example, 
k to the left (move r'). This move will induce a link with 
i with unit probability, but j will be added to the pseu- 
docluster with a probability depending on its energies of 
interaction with i and k. We shall require that the re- 
verse move be initiated by the seed particle that began 
the forward move. However, displacing i to the left ('re- 
verse' move r) will result in links with j and k only with 
probability 1 — g-^^^^" , where a € {j, fc} and Ae^ is the 
change in binding energy between i and a. To ensure 
that rightward and leftward moves of C occur with equal 
probability, given that the seed particle is i, we must sup- 
press the likelihood of the forward (rightward) move by a 
factor (1 — e~^^'^^'){l — e"^^'^"'"). For strong attractions 
(the regime in which we are interested) this factor is close 
to unity; for hard particles with no attractions it is zero, 
and we must reject the forward move [26]. We thus en- 
sure that the probability of a given particle 'pushing' or 
'pulUng' its host cluster is identical [57] • 

The acceptance probability for the virtual-move pro- 
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FIG. 4: Ensuring reversibility for cluster moves. We require 
that the likelihood of a given particle, say i, 'pushing' (f) and 
'pulling' (r) its host cluster is identical. To do so, every time 
we link two particles we record the likelihood that the link 
would have formed had the reverse move been proposed. In 
the reverse move, the linking particle translates in the oppo- 
site direction (and/or rotates with the opposite sense). We 
account for this probability difference via the overall accep- 
tance rate, ensuring superdetailed balance. The acceptance 
rate for such collective moves is high if the binding energy 
of particles is high (the regime of interest), and vanishes in 
the limit of vanishing interaction strength [26]. We require 
that the reverse move be initiated by the seed particle of the 
forward move, because moves having distinct seeds (compare 
f and r') cannot in general balance each other. 

cedure follows from Equations (|5| and ^ . We first eval- 
uate the link- failure factors qij. Consider a move from 
state PL to state v. The probability of not linking a parti- 
cle j (that is not a member of the chosen pseudocluster) 
to a particle i (a member of the pseudocluster) is, from 
Equation ([6]), 

q,,(^ ^ z.) = min {l.e^^^'-^^) (9) 

Here e|j^ denotes the pairwise energy of the bond ij in 

state a. Recall that the factor t['^^ is unity if in state 
/i particles i and j lie close enough to interact, and zero 
otherwise. 

The factor corresponding to ^ for the reverse move 
follows by interchanging the labels /i and u. Hence we 
have for external failed links 

where {ij)' denotes all shaded-white pairs except those 
which fall into two classes. The first class consists of 
those particle pairs that do not interact in state ^, but 
which move uphill in energy upon going from state p, to 
state u. Since we define noninteracting particle pairs (via 
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the interaction criterion I'j''') as those that in state /i lie 
too far apart to possess any energy of interaction, mov- 
ing uphill in energy means that these pairs end the move 
(state ly) with positive energy (we shall refer to particles 
having positive energy of interaction as 'overlapping' par- 
ticles) . The second class of particle pairs are those that 
move downhill in energy upon going from state fi to state 

but do not interact in state v. For the interaction cri- 
terion as defined this means those particle pairs that in 
state n overlap each other, but which finish their move 
(in state v) outside the interaction region. Note that 
the second class cannot exist if the potential is composed 
of a purely attractive piece plus a hard-core repulsion, 
because then no particle pairs start in an overlapping 
position. 

Our final acceptance rate is then 

w^acc(M ^i^\n) = 

e (Tie - nc)V{C) min J 1, e~''(''^'''"'-^'") 



(11) 



The link-failure factors external to the pseudocluster 
have canceled the ratio of Boltzmann bond weights, ex- 
cept those corresponding to particle-pairs {ij)n^o falling 
in either of the two classes defined above: class 1 contains 
those particle pairs that start (fi) in a noninteracting con- 
figuration and end (ly) in an overlapping one; class 2 con- 
sists of those pairs that start (/i) in an overlapping con- 
figuration and end (z^) in a noninteracting one. The sub- 
script 'n^o' stands for 'noninteracting<->overlapping'. 
Link- failure factors internal to the pseudocluster (the sec- 
ond product on the right-hand side of Equation (111) 
must be accounted for exphcitly. 

The third product on the right-hand side of Equa- 
tion (11 1 runs over one particular realization TZ of links. 



It can be thought of as a means of ensuring that the 
probability of a given particle 'pushing' a cluster is the 
same as its probability of 'pulling' the same cluster (or 
for rotations, ensuring an equal likelihood of effecting a 
rigid-body rotation in clockwise and anticlockwise direc- 
tions), and is required to enforce balance. This factor 
approaches unity for strongly-bound particles, indicat- 
ing that in this regime multi-particle displacements oc- 
cur with high probability. The factor becomes small for 
particles attracting weakly. 

Returning to our earlier example. Figure |3] we see that 
the only possible contribution to the first factor (of the 
second argument) of the acceptance criterion ( |11[ ) can 
come from bond jm (because in state /i this bond is 'non- 
interacting' according to our criterion). However, we see 
that in the proposed new state particles j and m do not 
overlap, and so their pair energy does not enter this first 
factor. The latter is therefore unity. Considering the fi- 



ijk possesses a large binding energy the proposed collec- 
tive displacement is likely to be accepted. In this case 
we accept state v according to our choice of its diffusion 
properties. 

One strength of the virtual-move scheme is that it al- 
lows one to identify groups of particles that move in con- 
cert, and to suppress their displacement in order to ap- 
proximate hydrodynamic damping. A body moving in 
an overdamped fashion through a fluid possesses trans- 
lational and rotational diffusion properties that depend 
on its size and shape. These properties cannot be fully 
captured by considering only individual or pairwise inter- 
actions between the cluster's constituent particles. Sim- 
ulating the solvent flow that mediates such many-body 
forces is extremely expensive computationally, and so col- 
lective hydrodynamic effects are often neglected when in- 
tegrating Brownian equations of motion. 

Stokes' law states that a cluster of effective hydrody- 
namic radius R (a measure of the greatest extent of the 
cluster perpendicular to the direction of motion) pos- 
sesses a translational diffusion constant Dt{R) = r/(3ii), 
and a rotational diffusion constant -Di(i?) = T/{8R'^), 



where F = k^T^-Kj^w) 



Vn 



10 ^ Pa s is the vis- 



cosity of water. We can respect this damping within 
the virtual-move algorithm by calculating for the chosen 
pseudocluster its effective hydrodynamic radius R, where 

(i?-i?o)'-(|(r.-re)xn|2). (12) 

i?o is the monomer radius. The average (■) = n^^ X]"^! 
runs over all of the nc particles i (having coordinates r^) 
comprising the pseudocluster (in a particular conflgura- 
tion) . The vector n is either the direction of translation 
or the axis of rotation, as appropriate; the vector Tc is for 
rotations the position of the centre of rotation, and for 
translations the centre of mass (r^) of the diffusing pseu- 
docluster. This size- and shape-dependent drag becomes 
increasingly important when the system in question is 
composed of very poly disperse or anisotropic aggregates. 

We enforce this damping by suppressing cluster dis- 
placements by a factor 2?-y(i?), where 7 S {t,r} for a 
translation or rotation as required. We set 'Dt{R) = 
Ro/R, and V,{R) = {Ro/Rf. 

We consider in Section |IV| a system of particles with 
purely attractive interactions and hard-core repulsions, 
in which case the second class of particle pairs in the 
set (ij) n^o does not exist: particles may not start in an 
overlapping position. In this case Equation (111 reduces 
to 



W^cc ifJ- ^ i^\n) = Q{n^~nc)V^{R) 



n 

X mm < 1, JJ^ e^*^'" 



q.jifi v) 



n 

n 



(13) 



nal two factors of Equation (111, we see that if the trimer Here the product n(ij)„^o i'™^ over pairs that do not 
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interact in state ^ (so that e^^^ = 0) but possess positive 
energy (overlap) in state u. Barring such overlaps this 
factor is unity; if such overlaps occur then this factor is 
zero, and we reject the move. For potentials permitting 
'soft' overlaps, Equation (111 does not reduce to Equa- 



tion (13 1, and overlaps are rejected probabilistically. 



The function Q is used in conjunction with the early 
rejection scheme (see Equation (|6|). This scheme effects 
a suppression of the generation rate of moves of clusters 
of size ric by a factor of 1 /nc , ensuring that all particles 
move with approximately equal frequency. 



In Section III we shall evolve a system of Lcnnard- 
Joncs discs using both VMMC and a simple Brownian 
dynamics protocol that neglects collective hydrodynamic 
effects. The latter enforces a translational diffusion con- 



stant for a cluster of size nc that scales as 



and a 



rotational diffusion constant that scales as the reciprocal 
of the cluster's moment of intertia. Because disc-disc in- 
teractions are isotropic, simply rotating a seed particle 
about an axis through its centre cannot induce a col- 
lective rotation. We therefore use as our basic move a 
combination of a translation and a rotation. To ensure 
that the resulting collective translational and rotational 
diffusion behaves as it would under Brownian dynamics, 
we generalize the procedure of this section to one that 
permits the use of distinct real and virtual moves. 

This completes our discussion of the the key result of 
this paper, the virtual-move Monte Carlo scheme. By 
forming links according to individual bond energies be- 
fore and after a proposed move, we displace particles 
collectively according to individual bond energy gradi- 
ents without calculating forces explicitly. We ensure that 
particle positions are updated with approximately equal 
frequency, and we damp the movement of multi-particle 
clusters in order to respect Stokes' law. By doing so, we 
can restore to the Monte Carlo procedure the collective 
diffusive motion suppressed by making sequential moves 
of particles in the face of strong, short-ranged interac- 
tions. In the following section we test this algorithm 
against Brownian dynamics simulations. In Section |IV| 
we apply VMMC to a model of biological self-assembly. 



III. AVOIDING UNPHYSICAL KINETIC 
TRAPS IN THE FACE OF STRONG 
INTERACTIONS 

Application to a schematic model of aggregation. Be- 
fore applying the virtual-move scheme to an example of 
self-assembly, we first verify that it evolves a system of 
particles according to an approximation of natural dy- 
namics. We consider a two-dimensional system of 324 
discs of radius a. Pairs of discs whose centres are sepa- 
rated by a distance r interact via a Lennard-Jones poten- 
tial modified to effect a range of attraction that is short 
compared to a particle's size: 



Here Tc = 2.5a and eb = 50 fceT" are respectively the 
range and strength of the interaction, and f = r — a de- 
notes a shifted distance. We have introduced a 'Lennard- 
Jones' function C{x) = A{x~^'^ — x~^). The second term 
in Equation ( [l4| ) 'shifts' the potential to zero at a cutoff 
distance of Vc', as a conseqence, the potential minimum is 
approximately —S^k^T (instead of —^Ok^T). Particles 
occupy about 10% of the box area. 

We shall evolve this system according to a simple 
Brownian dynamics protocol in which individual parti- 
cles experience a random force but not a random torque. 
Collective hydrodynamic effects are ignored. We update 
particle positions according to the set of equations 



(15) 



where i labels particles, the forces F; are derived from 
the potential (14), and 7 is a friction coefficient. The 



term 77^ is a Gaussian random force with correlations 



{'n^{tH{t')) = 2kBTjS{t - i')%l- 



(16) 



u(r) = ebO (re - r) [Cir/a) - Cif,/a)] 



(14) 



Such a dynamics effects local motion according to poten- 
tial energy gradients. It also promotes a buffeting-driven 
rigid-body translational diffusion of clusters of size nc 
scaling as n'^^ , together with a rigid-body cluster rota- 
tion scaling as 1^^ , the reciprocal of the cluster's mo- 
ment of inertia about the rotation axis (relative to that 
of a monomer) . Note that these results differ from those 
implied by Stokes' law. 

Our aim in this section is to use the virtual move 
scheme to mimic this dynamics. Because pair interac- 
tions are isotropic, collective rotations may not be initi- 
ated by rotating a seed particle about an axis through its 
centre. We therefore use as our basic move a combination 
of a translation and a rotation. To ensure that the emer- 
gent collective motion scales in the appropriate manner 
it is necessary to generalize the virtual move scheme de- 
scribed in the previous section to one in which we use 
distinct virtual and real moves. We shall use a virtual 
move consisting of the superposition of a translation and 
a rotation to select from our system a pseudocluster. We 
shall then make a real move of the pseudocluster which 
consists of a translation and rotation in general different 
from those of the virtual move. 

The acceptance rate for such a procedure is corre- 
spondingly different from Equation (|11[ ). To explain its 
form, we use as an illustration Figure [5| Here we consider 
particles i, j, and fc, all of which interact. Distances are 
exaggerated for clarity. From starting state /i the seed 
particle i is assigned a virtual move map consisting of 
the superposition of a leftward translation and a rota- 
tion about the z axis through the seed's centre (in our 
two-dimensional example we consider particles to move in 
the X — y plane) . Virtual translation magnitudes and ro- 
tation angles are drawn from uniform distributions with 
respective maxima A = O.llcr and A^t ~ 10°. We pro- 
pose links to all particles with which the seed interacts in 



9 



o % « o ^ 




(1) 



FIG. 5: Generalizing VMMC to allow distinct real and virtual 
moves. From starting state /x the seed particle i is assigned 
a virtual move defining a leftward translation and a rotation 
about the z axis through its centre (here particles are confined 
to the X — y plane). In this example the link ij forms but ik 
and jk do not; we therefore record the link- forming probabil- 
ity pij(/i — > if) and the link-failure probabilities q{i,j)k{t^ — > if) 
[frames (1) and if]. Frame if (see text) shows the intermedi- 
ate (virtual) state for the forward move. We then return the 
chosen pseudocluster to its original position (not shown) and 
execute our chosen real move. Our real move consists of the 
translational and rotational components of the virtual move 
with the rotation diminished by the square root of the clus- 
ter's moment of inertia about the relevant axis. Here this 
move yields proposed final state v. From v we record the 
probability of making link ij by executing the reverse virtual 
move, Pij{p ir), and the probabilities of failing to form 
links ik and jk [frames (2) and ir]. Frame ir shows the inter- 
mediate (virtual) state for the reverse move. The acceptance 
probability for this procedure is given by Equation 



state /i. In this example the link ij forms but ik and jk 
do not; we therefore record the link-forming probability 
Pij{lJL if) and the link-failure probabilities (?ifc(/i — > if) 
and qjk{^J' — ^ if) [frames (1) and if]. We refer to the 
state formed by application of the forward virtual move, 
starting from state ^, as the intermediate state for the 
forward move, if. With the linking scheme finished, we 
return the pseudocluster ij to its original position (not 
shown) and execute our chosen real move. 

Our real move also consists of a superposition of a 
translation and a rotation, but the latter is performed 
about the z axis through the centre of mass of the pseu- 
docluster. This allows precise control of rotational and 
translational degrees of freedom. The direction and mag- 
nitude of the translation are identical to that of the vir- 
tual translation. We shall account for the diminished 
diffusion constant of the cluster via the early-rejection 
scheme, discussed below. We obtain the rotation angle 
by scaling the magnitude of the virtual rotation by a fac- 
tor {Ic /nc)~^^'^ ■ It is also possible at this stage to scale 



cluster translations and rotations in order to approxi- 
mate a hydrodynamic damping. Here we instead employ 
a scaling designed to mimic Brownian dynamics. 

In our example. Figure [5] the real move produces pro- 
posed final state v. We then ensure reversibility by exe- 
cuting a reverse virtual move from state v, and recording 
the link-forming probability Pij{i^ —> ir) and the link- 
failure probabilities gife(i^ ir) and qjkiy ir)- The 
reverse virtual move consists of the forward virtual move 
with sign inverted (translations occur in the opposite di- 
rection; rotations are performed with opposite sense to 
the forward rotation, but about the same axis relative 
to the seed particle position). Frame ir shows the in- 
termediate (virtual) state for the reverse move. These 
factors of p and q ensure that the likelihood of passing 
from ^ to V and back again, via intermediate state if in 
the forward (/i — s- v) direction and intermediate state ir 
in the reverse direction (y —f /i) is such that balance is 
preserved. The procedure of balancing forward and re- 
verse rates for a particular realization of pseudocluster 
links and for specified intermediate virtual states could 
be termed 'superduperdetailed balance'. 

For the case of distinct real and virtual moves we have 
instead of Equation (ITT]) 



M^acc(Ai ^ v\n) = 
Q (uc — nc) min < 1 



r) -Q Pij{l^ ir) 



(17) 



The cutoff Uc is drawn from the distribution Q{nc) oc 
One factor of ensures that particles move with 
roughly equal frequency. The second factor accounts for 



the fundamental diffusion rate of i 



for a cluster of size 



^ . This procedure represents a coarse-graining of the 
dynamics (a cluster of size nc moves one distance unit ev- 

— 1/2 

ery nc sweeps, instead of n^ distance units per sweep), 
and offsets to a considerable degree the waste associated 
with building large pseudoclusters only to suppress their 
moves by a large factor. The Boltzmann bond weights 
include all interactions between the pseudocluster and 
its environment. The factors of q and p account for links 
unformed and formed, respectively, during both forward 
and reverse virtual constructions. Unformed links in- 
clude those internal to the pseudocluster. The choices 
of virtual and real moves determine the initial and pro- 
posed final states, ^ and and the intermediate states 
in either direction. 

This scheme allows us to perform cluster rotations and 
translations whose relative and absolute rates approxi- 
mate those induced by our Brownian dynamics proto- 
col. While seemingly complicated, the acceptance rate 



( 17 1 is straightforwardly evaluated during the virtual- 



move procedure. Furthermore, we find that for late-stage 
coarsening the acceptance rate is reasonable: for isolated 
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whole-cluster motion the Boltzmann bond weights and 
the products over the factors q reduce to unity; the link 
factors p are frequently close to unity, give the strong, 
short-ranged nature of the interaction (see Appendix B). 
Moves that result in the motion of single particles are 
unaffected by the rescaling. 

Starting from a high-temperature configuration we 
evolve the Lennard- Jones system according to either 
Brownian dynamics (BD), virtual- move Monte Carlo 
(VMMC) or single-particle Monte Carlo (SPMC) pro- 
tocols. We show in Figure |6] configurations as a function 
of time obtained under all three protocols. We compare 
Brownian dynamics and Monte Carlo timescales by re- 
porting times in units of t^, the characteristic time for 
a free monomer to diffuse a length equal to its diame- 
ter. We observe similar behaviour as a function of time 
for BD and VMMC protocols, with the strong, short- 
ranged interaction inducing clustering. Clusters diffuse 
and merge, leading to non-compact, kinetically frustrated 
structures |28| . However, single-particle moves (using the 
same displacement scale) suppress unphysically the diffu- 
sion properties of clusters, leading to much smaller struc- 
tures on equivalent timescales. 

For the system considered here the computational ef- 
fort required by VMMC is less than that required by 
the Brownian dynamics protocol. For the latter, a small 
integration time step is required to maintain numerical 
stability in the face of the strong, short-ranged inter- 
action; consequently, about 1.2 xlO^ integration steps 
are required to advance the system one unit. By 
contrast, within the VMMC framework we can employ 
much larger basic steps, limited only by the requirement 
that particles take many steps to traverse typical inter- 
particle distances [37] . For our chosen distribution of dis- 
placements (whose magnitudes are drawn uniformly from 
the interval [0,0.11(7]) we require only « 1000 sweeps 
to advance the system one td unit. For the trajecto- 
ries shown, the processor times required to advance the 
system (250,3500,6000) id units are « (3.5,84.5,152.0) 
hours for BD, and « (0.5,6.0,10.0) hours for VMMC 
(processor times reported to the nearest half-hour). 

We estimate that our chosen VMMC displacement 
magnitude is as large as one can employ for this system 
and density. The clusters formed under VMMC were 
(for times t < 100 ^d) slightly more 'ragged' than those 
formed under BD. Thereafter, we could not distinguish 
clusters formed by the two protocols (cluster morpholo- 
gies differ considerably between different trajectories un- 
der the same algorithm). Thus, while our implementa- 
tion of VMMC is not a perfect reproduction of BD, it 
represents a good approximation thereof. Our calcula- 
tions (Appendix B) indicate that SPMC would retain a 
reasonable acceptance rate (and therefore begin to ap- 
proximate Brownian dynamics) only when the displace- 
ment maximum is reduced to less than A — cr/100. This 
reduction would diminish the timescale corresponding to 
a single Monte Carlo sweep by a factor of at least ~ 100, 
and would lengthen simulation times by the same factor. 



The considerations of Appendix B indicate that the 
degree to which VMMC can effect collective motion at a 
controllable rate depends upon basic displacement scale 
(for both real and virtual moves), and must be assessed 
carefully for the model under study. We speculate that 
the utility of VMMC is greatest for models with short- 
ranged and anisotropic interactions, for models whose 
simulation requires a large basic step size (e.g. models 
defined on a lattice), and for systems whose components 
are present at low concentration. 

Application to hard particles. The ability to use dis- 
tinct virtual and real moves can be used to make collec- 
tive moves of hard particles without attractive interac- 
tions. A suitable algorithm is as follows. We dictate that 
forward and reverse moves are initiated by the same vir- 
tual move (translation or translation plus rotation) of a 
given seed particle. The real move consists of either the 
virtual move (with probability 1/2) or the virtual move 
with sign inverted (with probability 1/2). In the latter 
case the real translation occurs in the direction opposite 
the virtual translation, and the real rotation possesses 
the opposite sense of the virtual rotation. In this man- 
ner a given particle can both 'pull' and 'push' its host 
cluster, regardless of the strength of energetic interac- 
tions. As an example, the three-member cluster shown 
in Figure |4] may be moved both rightwards and leftwards 
via the virtual move depicted in frame (f). This scheme 
allows one to execute 'avalanche' moves of hard particles 
without regard for the geometry of the avalanche [55] . 



IV. APPLICATION OF VMMC TO 
SELF-ASSEMBLY 

Self-assembly is the process whereby interacting com- 
ponents organize spontaneously into thermodynamically 
stable patterns or aggregates. It is the means of forma- 
tion for many biological structures, including the protein 
packaging of viruses [22^ i23j to the lipid membranes en- 
closing cells. The self-assembly of nanometer-scale ob- 
jects constitutes an important branch of nanotechnology. 

In studying self-assembly, it is important to identify 
the nature of the inter-component interactions that can 
lead to stable structures with the required symmetry. 
However, a full understanding of how assembly is effected 
requires a characterization of its dynamics. Binding that 
depends on the precise alignment of neighbouring com- 
ponents might ensure thermodynamic stability, but the 
time required for two bodies to collide in this manner 
could be prohibitively long. Similarly, very strong con- 
tacts contribute to the stability of equilibrium structures 
but could also transiently stabilize malformed aggregates, 
thereby impeding equilibration. 

Brownian dynamics is perhaps the most natural kinetic 
model for several examples of biological self-assembly, 
given the relatively large sizes and small diffusivities of 
many proteins. However, individual components diffuse 
much more rapidly than do large-scale structures. A 
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FIG. 6: Configurations as a function of time for very attractive discs evolved according to Brownian dynamics (BD), virtual- 
move Monte Carlo (VMMC) and single-particle Monte Carlo (SPMC) protocols. Cluster sizes and shapes as a function of time 
are similar under BD and VMMC algorithms. VMMC and SPMC simulations were performed using the same distribution 
of basic displacements. The acceptance rate for sequential moves of single particles in this regime is low enough to suppress, 
unphysically, collective modes of relaxation. The computational efficiency of VMMC exceeds that of the BD protocol by more 
than an order of magnitude. 



faithful accounting of the faster of these two motions re- 
quires, in the face of strong, short-ranged interactions, a 
small integration time step in order to maintain numeri- 
cal stability. Long simulation times are therefore required 
in order to study the assembly of overall structures. 

Here we demonstrate that the virtual-move Monte 
Carlo algorithm can be used to evolve in a computa- 
tionally efficient manner a collection of self-assembling 
components possessing strong, short-ranged and angu- 
larly specific interactions. The chief advantage of VMMC 
over Brownian dynamics when applied to such models is 
that we may use with the Monte Carlo protocol a basic 
translation step that is not restricted by the width and 
depth of the pair potential well. With VMMC we face 
only the less stringent restriction that many steps must 
be taken to traverse typical inter-particle distances [37] . 

In the presence of ATP and magnesium ions the heat 
shock protein (Hsp60) from the organism Sulfolobus shi- 
batae self-assembles in two stages. First, monomers ('sub- 
units') of the protein assemble into 18-membered, nearly- 
spherical complexes ('units') of radius 17 nm. Next, 
units aggregate into extended structures often microns 
in scale. We refer to the 18-membered complexes as, in- 



terchangably, 'units' or 'chaperonins'. Experiments with 
the wild-type protein lead to two distinct types of chap- 
eronin super-structures under subtly different external 
conditions: two-dimensional sheets with a high degree of 
hexagonal order, and quasi one-dimensional strings. We 
focus here on assembly into sheets, the control of which 
provides a means of engineering organic templates with 
a high degree of order, potentially useful for electronics 
devices. Free-floating sheets are also formed by the self- 
assembly of inorganic nanocrystals [5^ . 

Computer models can help reveal the range of inter- 
unit attraction strength and specificity required to effect 
large-scale self-assembly [33l E]. The detailed interac- 
tions between chaperonins are not known. Patterns of 
polar, nonpolar and charged amino acid side chains ex- 
posed by individual chaperonins do not immediately sug- 
gest specific regions where two units would strongly bind, 
nor do experimental results over a range of ionic strength 
clarify the physical nature of the binding interaction. It is 
clear, however, that forces between sheet-forming chaper- 
onins favour equatorial contact [36] . We can exploit this 
information when constructing a simple model of chap- 
eronin self-assembly. 
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a parameter a) of the value a. The symmetrized cooper- 
ativity function Cq {^p) = Ca (V')+C-q {ip) rewards values 
of cos tp near ±a. The first angular factor in ( 18 1 encour- 



FIG. 7: Geometry for a schematic model of chaperonin self- 
assembly [161 1171 I18| . We use a short-ranged, anisotropic 
pair potential to mimic the tendency of chaperonins to bind 
equator-to-equator. We show chaperonin structure deter- 
mined by homology (image courtesy of Matthew B. Francis 
and Chad D. Paavola) together with the angles relevant for 
our chosen interaction. The angle between orientation vectors 
is 0; the angles between orientation vectors and the inter-unit 
separation vector are 9i and 82. 



We build such a model by coarse-graining over the mi- 
croscopic details of individual units (we consider units 
to be stable against dissociation into protein monomers) . 
The simplest such approximation is to regard units as 
spheres without surface detail [35j. We mimic the ef- 
fect of the microscopic unit-unit interactions by endowing 
spheres with a short-ranged, anisotropic pair potential 
designed to encourage mutual equatorial contact. 

Figure [7] illustrates the geometry of our chosen poten- 
tial. Each sphere has a polar orientation vector, shown 
as a double-headed arrow (we assume 'up-down' sym- 
metry). We also assume azimuthal symmetry, so that 
the pair potential does not change if spheres are rotated 
about their orientation vector (chaperonin units possess 
9-fold rotational symmetry about their polar axis). We 
choose the attraction to be strong if neighbouring orien- 
tation vectors are aligned (0 = 0, tt), and perpendicular 
to the inter-unit separation vector {61,2 = tt/2,3tt/2). 
Thus we allow only binding between complementary re- 
gions (equator-to-equator), and do not allow, for ex- 
ample, equator-to-pole binding. This reflects the near- 
perfect alignment observed between chaperonins in sheet- 
like assemblies, and also the notable absence of disor- 
dered aggregates 

Analytically, our chosen potential is 

e(r) = - Jeq e (i?„,ax - f) C,{^) Co{9i)Co{e2). (18) 

The step function O ensures that two spheres interact 
only when their surfaces are separated by a distance less 
than i?max, which we shall vary between i?o/4 and Ro/8, 
Rq being the chaperonin radius. Here r = r — 2R > is 
the distance between the surfaces of neighbouring chap- 
eronins. For f < we assume a hard-core repulsion. We 
parameterize the angular interaction via the 'cooperativ- 
ity' function Ca {ip) = e"'™*"^"") /'^ , which rewards an 
angle ■0 if its cosine is within some tolerance (specified by 



ages the alignment of neighbouring units; the second and 
third angular factors encourage mutual equatorial con- 
tact. We vary a between 0.2 and 0.4. 

Here we demonstrate that the virtual-move algorithm 
(incorporating hydrodynamic damping) described in Sec- 
tion |TT] can identify several distinct mechanisms of sheet 
assembly. These range from monomer addition to a single 
growing cluster (a nucleus), to the binding and merging 
of separate sheet-like structures. The latter mechanism, 
which turns out to be a plentiful source of kinetic 'traps', 
is strongly suppressed, using the same displacement scale, 
by conventional single-particle MC methods. We focus 
here on typical mechanisms of assembly for our model 
system, and leave a detailed study of the rates of ag- 
gregation as a function of attraction parameters, as well 
as the mechanisms underlying the formation of hybrid 
sheet- and string-like structures, for elsewhere. 

By assuming a chaperonin radius of i?o = 9 nm, we 
obtain from Stokes' law a translational diffusion con- 
stant for chaperonin monomers of Dt{Ro) w 5 x 10^^^ 
m^s""'^ (compare the self-diffusion coefficient of water, 
^'hsO ~ 5 X 10~^ m^s~-^). We draw displacements 
from a uniform distribution with magnitude 0.9i?o- Our 
monomer displacement timescale corresponds approxi- 
mately to 10^^ s. We fix the relative rates of transla- 
tion and rotation by imposing Stokes' law for monomers, 
namely 3 Dt{Ro) — 8Rq D,-{Ro). We ran simulations for 
^ 10^ to ^ 10^ MC sweeps, and so probe timescales of the 
order of seconds. In experiments [ini |T7l [T5] , large-scale 
assembly is observed on timescales of minutes to hours. 
Thus we expect our dynamic simulations to detect at 
least the onset of significant chaperonin self-assembly. 

We start from an initial state consisting of 1000 
monomers randomly dispersed and oriented in a three 
dimensional simulation box with periodic boundaries in 
each dimension. Units comprised about 0.8 % by vol- 
ume of the simulation box, equivalent to a protein con- 
centration of about 5 mg/ml (experimental concentra- 
tions of protein range between 1 and 5 mg/ml |M])- We 
then evolve the system according to the virtual-move 
scheme with a hydrodynamic damping. Times are quoted 
in virtual-move Monte Carlo sweeps, with one VMMC 
sweep corresponding to the (uniform) choice of 1000 seed 
particles. 

In Figure [8] we present a 'kinetic phase diagram' of 
chaperonin assembly, obtained from single trajectories. 
We vary interaction strength J^q and inverse angular 
specificity a. For interaction ranges i?rnax — Ro/4: (top 
panel) and i?niax = ^0/8 (bottom panel) we indicate 
where we find 'no assembly' (open squares), 'good as- 
sembly' (circles) and 'bad assembly' (closed squares). We 
denote by tend the largest time accessed at each param- 
eter point (the 'end' of the trajectory). We conclude 
that 'no assembly' has taken place if the largest cluster 
at time tend possesses fewer than 15 chaperonin units. 
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FIG. 8: Kinetic phase diagram for our schematic chaperonin 
system, evolved using the VMMC algorithm with hydrody- 
namic damping. We indicate regions of 'no assembly' (open 
squares), 'good assembly' (circles) and 'bad assembly' (closed 
squares); see text. A good assembly-bad assembly pair indi- 
cates that intermediate building blocks are well-formed, but 
subsequent multi-particle collisions induce kinetic frustration. 
Configurations corresponding to parameter sets a, b and c are 
shown in Figures [9lto |ll[ 




e 




o 






FIG. 9: Configurations obtained after 1.2 x 10® (top) and 
8 X 10*^ (bottom) VMMC sweeps for our schematic chap- 
eronin system of spheres with sticky equators of strength 
Jcq = 6.5fcBr (parameter set a of Figure [8|. At this attrac- 
tion strength nucleation is sufficiently rare that self-assembly 
proceeds by the binding of monomers to a single sheet. The 
resulting structure is relatively well-formed. 



For those systems for which this is not true, we conclude 
that assembly is 'good' if the constituent monomers of 
the largest cluster in the system possess on average more 
than 4.75 'bonds'. We define a particle's bond number 
as its energy divided by the equatorial coupling, Joq. We 
conclude that assembly is 'bad' if, at tend, the system's 
largest cluster possesses more than than 15 members, but 
fewer than 4.75 bonds per member. 

For those parameter sets exhibiting bad assembly, we 
identify the largest number of bonds possessed by mem- 
bers of the largest cluster at any time along the trajec- 
tory. If this number is greater than 4.75, we plot a 'good 
assembly-bad assembly' symbol pair. This indicates that 
while the intermediate building blocks may at some time 
be well-formed, collisions between these multi-particle 



structures eventually give rise to an aggregate that is 
ill-formed. 

A note on timescales: ideally, we would present in a ki- 
netic phase diagram data obtained at fixed Monte Carlo 
time. Here that is not feasible, because of the broad 
distribution of relaxation times observed in our chaper- 
onin system. For example, for parameter set (-Rmax = 
i?o/4, Joq = SfceT, cr = 0.4) we observe the aggregation 
of all particles in ^ 10® VMMC sweeps. After an equal 
time, parameter set (-Rmax = -Ro/8, Jcq = 6.5A:bT, cr = 
0.3) exhibits no appreciable assembly (nucleation is very 
slow). At longer times, however, the latter system is well- 
assembled. Consequently, we show data at fixed proces- 
sor times, such that we judged assembly for all parameter 
sets to be sufficiently far advanced that the ultimate fate 
of each system is clear. This criterion is clearly arbi- 
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FIG. 10: Configurations obtained using the VMMC algo- 
rithm apphed to the chaperonin system with equatorial inter- 
action strength Joq = Ik^T (parameter set b of Figure [8|. 
A modest increase in attraction strength promotes nucleation 
to a considerable degree. Assembly proceeds both by the ad- 
dition of monomers to single sheets (top panel, ~ 0.4 x lO'' 
sweeps), and via collisions of multi- particle sheets (bottom 
panel, ~ 6.3 x 10^ sweeps). Sheets often collide awkwardly, 
producing ill-formed structures. Here relaxation of these 
structures is sufficiently rapid that assembly is still 'good'. 



trary. However, we believe that the mechanisms of as- 
sembly revealed in this manner are qualitatively robust 
to variations in the criteria used to draw the kinetic phase 
diagrams. 

The general trend of assembly revealed in Figure |8] in- 
dicates that regions of good assembly occupy a relatively 
small region of parameter space. This region is defined 
by a balance between collision rates (controlled largely by 
density, attraction range and the specificity parameter ct) 
and relaxation rate (controlled largely by Joq), such that 
assembled structures form rapidly enough to be observed 
on the timescales simulated, but not so rapidly that they 
malform. If this optimal ratio ('good assembly', circles) is 
disturbed, we observe over-rapid growth leading to mal- 





FIG. 11: Configuration obtained via VMMC for our 
schematic chaperonin system with equatorial interaction 
strength Jcq = %k^T (parameter set c of Figure [8|. Nu- 
cleation is so rapid that multi-particle binding events occur 
frequently. The resulting structures fail to relax before en- 
countering other such stuctures, leading to aggregates trapped 
far from equilibrium. Top: excluded-volume view. Bottom: 
bond view. 



formed structures ('bad assembly', closed squares), or 
growth too slow to be observed on the timescales sim- 
ulated ('no assembly', open squares). For one parameter 
set we see initially good assembly, where monomers bind 
into small, well-formed sheets, followed by bad assem- 
bly induced by sheets colliding awkwardly and produc- 
ing malformed structures. The rate of sheet collision is 
set by Stokes' law, which is respected by the VMMC al- 
gorithm. Reducing the attraction range (upper panel to 
lower panel) has the effect of reducing collision frequency. 
We observe an offsetting and a narrowing of the region 
of good assembly. 

It is illuminating to examine the assembly mechanisms 
observed as we vary only the strength of the equato- 
rial coupling, exemplified by parameter sets a, b and 
c in the top panel of Figure [s] In Figure |9] (parameter 
set a of Figure [s]) we show a late-time configuration for 
units possessing equatorial coupling Jcq = %k^T . As- 
sembly proceeds, following a rare nucleation event, via 
the binding of monomers to a single sheet-like nucleus. 
The resulting structure is well-formed and defect-free. 
The crossover from non-assembly to assembly is rather 
sharp: at the concentrations used, we observed no as- 
sembly within our simulation time at equatorial attrac- 
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FIG. 12: Configurations obtained via a single-particle Monte 
Carlo algorithm for our schematic chaperonin system with 
Rmux = Ro/4:, (J — 0.3 and (top to bottom) equatorial inter- 
actions Jcq = 6.5, 7 and 8 ksT. For the displacement distribu- 
tion employed here (uniform, with maximum displacement 0.9 
-Ro) single-particle moves strongly suppress collective modes 
of motion, and very few multi-particle collisions take place. 
For all attraction stengths structures are well-formed, even 
though VMMC results indicate that for the two larger values 
of Jcq a collective dynamics results in the formation of non- 
optimal aggregates via multi-particle collisions (see Figures |9] 
to [11]. 



tion strength Jcq = 5.5 k^T. 

The change in assembly mechanism caused by increas- 
ing the unit-unit interaction strength is dramatic. At 
a shghtly larger interaction strength, Jcq = 7 k^T (Fig- 
ure 10 parameter set b of Figure [s]), nucleation proves 
more rapid. Assembly proceeds via the organization 



of monomers into multiple sheets. These sheets diffuse 
according to Stokes' law, and collide with each other. 
Multi-particle collisions are often awkward, providing an 
ill-formed template to which monomers bind. In the 
example shown, however, relaxation is sufficiently rapid 
that large structures can relax via large-lengthscale fluc- 
tuations: assembly is still 'good'. 

At still higher attraction strengths, such as Jcq — 



SksT (Figure 11 parameter set c of Figure |8|, aggre- 
gation is so rapid that nuclei do not have time to relax 
into low energy sheet-like structures before they bind to 
other such ill- formed aggregates. In this regime assembly 
is frustrated kinetically. 

These results indicate that the assembly mechanism 
for our schematic chaperonin model changes consider- 
ably with the unit-unit attraction strength. Because 
these interactions are angularly specific, particles must 
collide equator-to-equator in order to bind. For insuf- 
ficiently strong equatorial attractions, random collisions 
between monomers do not result in stable intermediates, 
and no assembly is seen. For sufficiently strong attrac- 
tions, assembly proceeds via the binding of monomers 
to a single sheet-like nucleus. Aggregates grown in this 
way are typically we 11- formed and defect-free. Increasing 
unit-unit couplings beyond this point slows equilibration. 
Nucleation is promoted, and collisions between multiple 
sheet-like nuclei, which occur with a frequency governed 
by Stokes' law, usually result in awkwardly-bound struc- 
tures that relax only slowly. There exists a narrow regime 
of parameter space in which structures formed in this way 
can relax, via collective fluctuations, into an approxima- 
tion of a well-assembled sheet. However, at very large 
attraction strengths nucleation is so rapid that the nu- 
clei themselves are ill formed, leading to disordered ag- 
gregates. 

For the basic displacement scale considered here, 
single-particle Monte Carlo techniques suppress unphys- 
ically the diffusion of multi-particle structures. When 
interactions are such that appreciable clustering of parti- 
cles develops, the dynamics of a system evolved accord- 
ing to single-particle Monte Carlo protocols does not sat- 
isfy Stokes' law. As a result, the source of kinetic traps 
whereby clusters collide and bind awkwardly is strongly 
suppressed. In Figure [12] we show configurations gen- 
erated by single-particle translations and rotations ap- 
plied to systems corresponding to the interaction pa- 
rameters of Figures [9] [lO] and 11 (parameter sets a, b 
and c of Figure [8|. For all three attraction strengths 
single-particle moves generate we 11- formed, isolated clus- 
ters. For an attraction strength of Jcq = 6.5 k^T the 
VMMC results indicate that the ratio of the rates of clus- 
ter growth and diffusion is such that single-particle addi- 
tion to a single nucleus is the dominant assembly mecha- 
nism; single-particle moves naturally capture this dynam- 
ics. However, for the larger attraction strengths a dy- 
namical protocol satisfying Stokes' law generates multi- 
particle collisions, inducing a degree of kinetic frustration 
that increases with attraction strength. Single-particle 
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moves fail to identify this mechanism. Note that if inter- 
component interactions are such that assembly must pro- 
ceed via the interaction of multi-particle structures [25] . 
single-particle moves would encounter an unphysical ki- 
netic trap. 



V. CONCLUSIONS 

We have presented a virtual-move Monte Carlo clus- 
ter algorithm designed to permit the collective relaxation 
of particles possessing attractions of arbitrary strength, 
range and geometry, an important example being self- 
assembling particles endowed with strong, short-ranged 
and angularly specific ('patchy') attractions. By calcu- 
lating pair energies before and after notional ('virtual') 
moves, we deduce whether particles move individually or 
in concert. Using an 'early rejection' scheme designed 
to suppress moves of clusters by a factor inversely pro- 
portional to the cluster size, we ensure that all particles 
move with approximately equal frequency. We also en- 
sure that Stokes' law is satisfied for the explicit collective 
motion of aggregates of arbitrary size and shape. 

Our scheme approximates the simultaneous updates 
of particle positions characteristic of molecular dynam- 
ics simulations. Its advantage over the latter lies in the 
fact that by computing energies before and after vir- 
tual moves, one bypasses the need to compute forces and 
torques explicitly. In addition, one can propose particle 
displacements that are not limited by the width of the 
potential well, regardless of the depth of that well. The 
time savings thus accrued can be considerable (see Sec- 
and may be estimated generically as follows. 



nil 



tion 

Consider a system of particles evolving over some time 
r. We assume that over this period there exist un pair- 
wise interactions. The computational effort required by 
a Brownian dynamics simulation scales as 



Cbd ~ 



(19) 



where At is the integration time step and Cp is the cost 
of evaluating forces and torques for each particle pair. 
The integration time step must be such that a typical 
particle displacement effects a change in pairwise energy 
not more than k^T. We therefore have 



At 



1 / ecr \ ^ 
2D VJ) ' 



(20) 



where D is the particle diffusion constant, e is the width 
of the potential well in units of the particle radius cr, and 
J is the depth of the potential well in units of k^T (we 
assume for simplicity a triangular potential well instead 
of a square well) . 

The cost of a virtual-move Monte Carlo sweep scales 

as 



-VMMC 



MMC 



(21) 



where AiyMMC is the basic VMMC timescale, Ce is the 
cost of evaluating one pairwise energy, and a = 0(1) is 
a parameter accounting for the fact that we may require 
more than one energy evaluation per particle pair (e.g. 
when executing a 'reverse' virtual move in order to en- 
force superdetailed balance). The VMMC timescale (de- 
rived from the basic displacement scale) is not governed 
by the width and depth of the potential well, as would 
be the case for a single-particle Monte Carlo algorithm, 
but is determined by the less stringent requirement that 
many steps are taken over typical inter-particle distances. 
We set the typical displacement to a fraction / of the par- 
ticle radius. We estimate that 



At 



VMMC 



(Pacc/g) 

2D 



(22) 



where Pacc is the acceptance rate for particle displace- 
ments. Comparing Equations (19) and (21 1 reveals that 



c 



VMMC 



q;Ce V ^ 



(23) 



The savings associated with VMMC become more pro- 
nounced the stronger and shorter-ranged is the poten- 
tial. For the chaperonin model studied here we have 
/ - 1, J - 10, e - 1/5 and Pacc ~ 0.1 to 1. The ra- 
tio C-p / {oiC-Ej is of order unity (a « 2, but evaluating 
forces and torques is more costly than evaluating ener- 
gies). We obtain therefore Cbd/Cvmmc ^ 10 to 10'^. In 
this estimate we neglect some of the 'overhead' of the 
virtual-move procedure, associated with suppressing the 
generation rate for moves of large clusters (in order to 
perform updates of particles with approximately equal 
frequency). Nonetheless, for attractive interactions of 
short range this estimate indicates that we might expect 
considerable time savings using VMMC as opposed to 
Brownian dynamics. 

Further comparisons between the method presented 
here and Newtonian simulations must be performed be- 
fore the generic dynamical fidelity of the former can be 
determined. To this end, a study is underway \6V in 
which we compare Brownian dynamics with the VMMC 
algorithm, where each is used to effect the assembly of 
idealized protein capsomers into icosahedral virus cap- 
sids [23]. 

We expect that the virtual-move algorithm can be used 
to study the phase behaviour and aggregation mecha- 
nisms of a variety of self- assembling systems, both on- 
and off-lattice. 
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VI. APPENDIX A: USING COLLECTIVE 
MOVES TO FACILITATE RELAXATION 

Monte Carlo cluster algorithms [121 (TJ are used to 
evolve strongly-attractive particles in order to avoid the 
suppression of diffusive modes of motion that plague 
single-particle protocols. Clusters are identified and 
moved on the basis of properties of particles in the cur- 
rent state of the system, such as energy or degree of prox- 
imity. While leading to efficient diffusion of clusters, the 
internal relaxation of structures evolved in this manner is 
often under-represented: particles in close proximity or 
interacting strongly are liable to be moved collectively, 
and therefore will not re-arrange relative to each other. 
We demonstrate here that a straightforward modification 
of these algorithms can be used to efficiently relax, in a 
collective manner, strongly-interacting particles. How- 
ever, although equilibration can be facilitated, particle 
motion does not proceed solely according to local poten- 
tial energy gradients; the algorithm discussed here should 
therefore be regarded only as a scheme for sampling equi- 
librium ensembles. We discuss this point in more detail 
in Appendix B. The algorithm we describe is similar to 
that proposed by Troisi, Wong and Ratner [55] . 

We choose as the linking probability 



Pij ^ v) = max 



0,1 



exp 



(/3f"f( 



(24) 



Here (3i is a free parameter that functions as a fictitious 
reciprocal temperature. The term itf is a fictitious po- 
tential, and can be chosen for convenience. The simplest 
choice is to set the fictitious potential equal to the true 
potential. For some applications a convenient choice is 
Mt(e) = e0(ro — r): Uf(e) = e within some cutoff dis- 
tance tq, and zero otherwise. This renders potentials of 
arbitrary range amenable to the iterative linking scheme 
described in Section |ll] Note that Equation (24 1 depends 
only on the interaction energy between i and j in the 
initial state /i.This is a straightforward generalization of 
the form chosen by Swendsen and Wang [12] , and reduces 
to an off-lattice version of the Swendsen- Wang algorithm 
in the limit that the fictitious potential and temperature 
are set equal to their 'true' counterparts. 

Equation ^ gives as the acceptance rate for the move 

W^,,{li -^v)= min(l,cxp(-/3£;,,^ + /3f?7,,^)) . (25) 

Here E^^^ = E,y — and U„^^ = U^ — Uf_i, where Ea de- 
notes the interfacial energy between pseudocluster C and 



while C/q < denotes the attractive part of the fic- 
titious energy between C and its environment, Ua = 
J2x (0, Wf (cij)). We have for simplicity set the dif- 
fusion term 2? = 1 and taken the cutoff ric oo. 

In the case where the inter-particle potential is attrac- 
tive, eij < 0, and we choose a fictitious potential ut(ey) 
equal to the true potential eij, the procedure we have 
described is particularly straightforward to implement. 
Particles are linked according to a simple probability, 
Pij = 1 — e~''fl'^'jl, and the ratio of acceptance rates is a 
simple function of the change in energy resulting from the 
move. For the move fj, —i- v, the acceptance probability 



( 25 1 reduces to 



W^acc(M ^ J^) = min l^e 



.(A-/3)(i5. 



(26) 



For infinite fictitious temperature, /3t = 0, the likelihood 
of forming links between the seed i and any other parti- 
cle is zero, and so the algorithm executes single-particle 
moves with acceptance probability min (l, e"''^^''"^^'^). 
For a fictitious temperature equal to the true tempera- 
ture, /3f = /?, and for attractive interactions, the accep- 
tance probability is unity. 

By drawing a fictitious reciprocal temperature in the 
range (3f € [0,/3] we can interpolate between single- 
particle moves and rejection- free cluster moves. When 
applied to a tightly-bound aggregate of particles, the 
rejection-free scheme docs not allow the construction of 
pseudoclusters smaller than the aggregate. Thus the ag- 
gregate cannot relax through moves in concert of its con- 
stituents. In other words, the generation rate for many 
transitions involving the collective motion of tightly- 
bound particles is close to zero. By contrast, varying 
(3{ allows one to increase markedly the generation rates 
of these moves, at the affordable cost of reducing slightly 
their acceptance rates. As a result, one can choose from 
the aggregate pseudoclusters of arbitrary size, and thus 
propose collective internal relaxations in the presence of 
arbitrarily strong interactions. We refer to this procedure 
as cluster 'cleaving'. 



its environment in a given microstate, Ea 



E2 



In Figure 13 we demonstrate the advantages of the 
cleaving algorithm over the rejection-free cluster algo- 
rithm. We show example configurations from a system 
100 hard discs of diameter a in two dimensions, endowed 
with an attractive piecewise-linear pair potential of range 
a. With the inter-particle separation denoted by r, par- 
ticles experience a hard-core repulsion for r < a. The 
potential (shown in Figure 13 1 increases linearly from 
its minimum, eo = — 45fcBT', at r = a to — 15 k^T at 
r = 3a/4. From r = 3a/4 the potential increases linearly 
to zero at r = 2a. Thereafter, it is zero. This system 
possesses a thermodynamically stable ground state cor- 
responding to an hexagonal close-packed sheet. 

We find that relaxation is frustrated and that jammed 
structures persist for the single-particle moves in combi- 
nation with rejection- free cluster moves. These jammed 
structures owe their existence to the following mecha- 
nism. First, potential energy gradients encourage neigh- 
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-45 fceT < e < -44.5 A:bT. If e is outside this range, the 
fictitious potential returns zero. This ensures that only 
those particle that are optimally bound undergo signif- 
icant collective translations or rotations. The resulting 
dynamics is unphysical (see Appendix B) but leads to 
efficient relaxation of the system. Structures formed un- 
der the cleaving protocol are typically more compact that 
those generated by dynamical algorithms. 



VII. APPENDIX B: MONTE CARLO 
DYNAMICS 

A. Single-particle moves 

A Monte Carlo simulation consisting of a sequence of 
moves of individual particles can approximate natural dy- 
namics for systems in which relaxation is not dominated 
by diffusive modes of motion. As a simple example, con- 
sider a single particle in one dimension in a potential 
U{r). For sufficiently small basic displacement scale A, 
a Metropolis Monte Carlo trajectory is equivalent to the 
behaviour described by a diffusive Fokker-Planck equa- 
tion. 

Following Refs. [5J [3] we can demonstrate this equiva- 
lence for a particle at position r in a potential U(r). The 
master equation for the particle's motion is 



FIG. 13: Configurations from representative trajectories of 
our test system (see text) at times ti < t2 < t^. Left pan- 
els: single-particle moves plus rejection- free cluster algorithm, 
leading to a configuration trapped far from equilibrium. Right 
panels: single-particle moves plus cleaving algorithm. The 
latter circumvents the kinetic traps associated with rarely 
testing strong inter-cluster bonds. 



bouring particles to aggregate into dense clumps sepa- 
rated by voids of lower density, a relaxation that can 
be accomplished readily by single-particle moves. There- 
after, because of the strength of the inter-particle con- 
tacts, the rejection- free cluster algorithm fails to propose 
moves of one clump relative to another. Since the struc- 
ture does not readily relax in a single-particle fashion - 
the steeper gradient of the potential near its minimum 
disfavours the 'jumping' of single particles across a gap 
- the system finds itself in a kinetic trap. 

The cleaving algorithm can circumvent this trap. For 
sufficiently large T{ all contacts, regardless of strength, 
are tested. Here we draw the fictitious reciprocal tem- 
perature from a uniform distribution between and (3, 
P{f3i) = P^^Q{(3 — Pi), and we use a fictitious potential 



Uf(e) = 



(e < eo + S,) 
(e > eo + S,) 



(27) 



Here — 0.1 fee T is a cutoff energy. The fictitious 



potential ( 27 1 returns the true energy e if e is within 
a tolerance of the potential minimum, eq, i.e. if 



dtP{r;t) 



P{r';t)W{r' r) 
P{r;t)W{r^r'). (28) 



Here J^, = J^^^ dr' , P{r; t) is the probability of finding 
the particle at position r at time t, and W{r — > r') is the 
rate for moving from position r to position r' . For the 
Metropolis acceptance rate the latter reads 



The parameter Wq is 



r') = 14^0 Wgcn('^r) min (l, c~''[^(''''-^H 

(29) 

a reference frequency, and 
Wgen(<5'') = (2A)~^ the a priori probability for choos- 
ing from a uniform distribution a step 5r = r' — r va the 
range [—A, A]. The 'min' term encodes the Metropolis 
acceptance probabiHty. 



For sufficiently small A one can expand Equation ( 28 1 
in powers of f = r' — r. This is most easily done by 
changing notation from W{r' r) to W{r + f\ —f). The 
latter symbol means the rate for passing from configu- 
ration r -\- f {= r') to configuration r + f — f = r, and 
can be expanded in its first argument: W{r + f; — f) = 
W{r;—f) + fdrW{r;—f) -I- • • • . Likewise, P{r',t) — 
P{r, t) + fdrP{r, t) + - ■ ■ . To second order in dr, Equation 
(p8l reads 



1 



dtP{r;t) « dr {{r)P{r,t)) + -df. {{f^)P{r,t)) , (30) 
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with 



The master equation for the coordinate r is 



{f'') EE / dff* VK(r;-f). 

-A 



(31) 



The derivative-free term J^, P{r,t)W{r;'r) — 
Ir' ^ (^; t )W{r; —r) vanishes by symmetry. Equa- 
(|30| is a Fokker-Planck equation with drift velocity 



tion 



and diffusion constant D = 



Equation ( 29 1 we have to first order in f 



W{r;~f) 



Wo 
2A 



min(l, 1 + /3f C/'(r)) , 



(32) 



We have assumed that f|C/'(r)| ^ 1. 

The integrals (31 1 can be evaluated in a piecewise fash- 



ion, with the minimum function in Equation (32) return- 
ing its first argument when f > 0, and its second ar- 
gument when f < 0. To second order in A one can 
calculate that the drift velocity is proportional to the 
potential gradient, v = — (/3/6)A^VFo U'{r), and that the 
diffusion constant is D = A^M^o/3. This corresponds to 
a Langevin dynamics satisfying an Einstein (fluctuation- 
dissipation) relation —v/D — {(3/2)U'{r). Note that 
normally v is independent of temperature T (not cx P 
as here), and D oc T (not independent of T). This can 
be arranged [6J by making the fundamental attempt fre- 
quency Wo proportional to T. 

Thus Monte Carlo moves of single particles in a poten- 
tial take place in a dynamically realistic way, provided 
the basic step size is such that large changes in energy 
are not encountered. If large changes in energy are en- 
countered, drift and diffusion cease to be have the forms 
derived: the first correction to the diffusion term, for ex- 
ample, is oc WoA^U'{r). When faced with this problem 
one must either make the basic step size very small, in 
which case computational costs can be prohibitive, or re- 
cover diffusion by means of explicit collective moves. The 
algorithms we have introduced address this issue. In the 
remainder of this Appendix we shall show that the clus- 
ter cleaving algorithm (introduced in Appendix A) , which 
forms pseudoclusters by linking particle pairs according 
to their energies, corresponds to an unphysical dynam- 
ics (a dynamics in which the drift velocity is not simply 
proportional to the potential gradient) , and so should be 
regarded only as a scheme for sampling equilibrium en- 
sembles. By contrast, the VMMC algorithm, which links 
particle pairs in a manner consistent with their poten- 
tial energy gradients, corresponds to an approximation 
of realistic dynamics. 



B. Cluster cleaving algorithm 

Here we consider the dynamics of the cleaving algo- 
rithm described in Appendix A. We set the fictitious po- 
tential U{ equal to the true potential e. Let us consider 
the separation r = Xj — Xi between two otherwise iso- 
lated particles i and j, which interact via an attractive 
pair potential ^ij{r) < 0. 



dtP{r;t) = I P{r';t)W,{r' ^r;pf) 

P{r;t)W,{r^r';Pt). (33) 



Using Here /- = / d/Sf Q{l3{) J^^ df, where f ee r' - r; P(r; t) is 



the probability of observing a bond separation r at time 
t; and Wc{r r';f3{) is the rate at which the cleaving 
move changes the bond separation from r to r' . Recall 
that Q{f3f) is the distribution from which we draw the 
fictitious reciprocal potential. 

The rate Wc (the subscript 'c' stands for 'cleaving') at 
which the separation r changes follows straightforwardly 
from Equations ( 24 1 and ( 26 1 . We assume that either 



particle i or particle j is chosen as a seed, and displaced 
by the vector f. Then the bond separation r changes 
if 1) no link is proposed between i and j, and 2) the 
Monte Carlo acceptance probability is satisfied. Condi- 
tion 1) occurs according to Equation (24) with probabil- 
ity min [1, exp (/3feij(r))] = exp (/3f ejj(r)j. Criterion 2) is 
satisfied with a likelihood equal to the right-hand side of 
Equation (26 1. Hence 



Wc(r ^ r') = (2A)-^ Wo exp [/3fe,j(r)] 
xmin(l,e-('3-ft)N('-')--('-)]) 



(34) 



Once again. Wo is a reference frequency, and again we 
assume that we can expand the master equation ( 33 1 to 



second order in the small displacement A. We obtain the 
Fokker-Planck equation 

dtP{r; t) « -dr {v,sP{r; t)) + ^9, {DdrP{r; t)) , (35) 



with effective drift velocity 



The 'bare' drift velocity v is 
A2 - 
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W ■ip-Pf)dre,jir), 



(36) 



(37) 



(38) 



(39) 



which contains an integral over f3{ (this integral acts on 
any /3t-dependent factors to its right). By virtue of the 
position-dependence of ( 39 1 , the effective drift velocity is 



and the diffusion constant is 

A2 . 
D^-W. 

Here we have defined the position-dependent rate 
W=Wo f d/3fQ(/3f)e*^-W, 



not simply proportional to the negative of the potential 
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gradient. Instead, it is biased by a term depending on 
the exponential of the local bond energy: 



off = -^1^0/3 [^r£^Jir)] J dP{Q{Pt)e 



ftei,-(r) 



(40) 



This bias is not consistent with physical dynamics. The 
cleaving algorithm evolves the system according only to a 
rough approximation of true Langevin dynamics: a par- 
ticle's drift velocity is not simply proportional to the po- 
tential gradient it experiences, nor is the diffusion con- 
stant (38 1 independent of position. Only in the (unde- 
sirable)Timit of single-particle moves, Q(/3f) = S{Pf), is 
local physical dynamics restored. Note also that in the 
convenional rejection-free limit, Q{(3t) = S{f3{ — (3), drift 
and diffusion are both suppressed by a factor e^'^'^^'"', 
demonstrating that in the limit of strong attractions lit- 
tle motion is possible. 

The root of these difficulties is the fact that pseudoclus- 
ters are built according to pair energies, and not energy 
gradients. However, drift and diffusion still satisfy the 
required Einstein relation for evolution to equilibrium. 



D 



f3 



dr£ijir). 



(41) 



This condition is equivalent to the master equation ( 33 ) 
satisfying balance. 



C. VMMC algorithm 

Next we turn to VMMC algorithm described in Sec- 
tion |lT] With similar notation to before, the master equa- 
tion for the separation r between two otherwise isolated 
particles i and j is 

dtP{r-t) = / P{r';t)W,{r' ^r) 



(42) 



where J. = dr. The rate (the subscript stands 
for 'virtual') at which the separation r changes follows 
from Equation We assume that the virtual displace- 
ment of particles i relative to j results in a change in the 
bond separation r by f. This change is accepted if 1) no 
link is proposed between i and j, and 2) the Monte Carlo 
acceptance probability is satisfied. Condition 1) occurs 
according to Equation ^ with probability 

q,,{r -> r') = min [1, exp (/3e„-(r) - (3e,j{r'))] . (43) 

If condition 1) is satisfied then condition 2) is automati- 
cally true provided that no particle overlaps occur. Thus 



Wv(r 



Wo 
2A 



(44) 



In a regime in which the energy change induced by the 
basic displacement is small, we may expand Equation 



(42 1 in powers of A. This procedure yields a Fokker- 



Planck equation with a physically realistic drift velocity 
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and a position-independent diffusion parameter 
D = —Wo. 



(45) 



(46) 



The rates of drift and diffusion of the relative coordinate 
r per Monte Carlo sweep are twice the values given by 
Equations (45 1 and (46 1, since both i and j are selected 



as 'seed' on average once per sweep. These results agree 
with those obtained by assuming that i and j are subject 
to a Brownian motion described by the equations 



7ii = -d^.e{r) + rji; 
1^3 = -dx,e{r) + rij, 



(47) 



provided that the friction coefficient 7 is related to 
the Monte Carlo attempt frequency Wo by A^VFqT = 
6kBT. The Gaussian white noise satisfies {rji{t)rij{t')) = 
2kBTj6,jS{t-t'). 

We next consider the motion of the centre of mass R = 
^{xi + Xj) of the dimer ij. The centre of mass position 
may be changed if 1) no link is formed between i and j, 
in which case one particle moves relative to the other, or 
2) if a link is formed between the particles, in which case 
i and j may move collectively. The master equation for 
R is consequently 



dtPiR) 



Wo 



A 

fWo 



P{R)q,,{r ^ r') 



') min 1 



P 



^ J^^ P{R)p..,{r ^ r') min (^1, ^) . (48) 



The first two lines of Equation ( 48 1 arise from relative 
moves of i and j, while the second two lines describe 
moves of i and j in concert. The frequency parameter 
/ permits an adjustment of the rate of collective moves 
relative to those of single particles (it is controlled by 
the early-rejection scheme described in Section |ll|, and 
the scale factor s quantifies the chosen scaling of the col- 
lective displacement relative to that of a monomer. The 
'min' function enforces the requirement that forward and 
reverse collective moves must be initiated by a given seed 
particle moving in the forward and reverse directions; we 
have defined p = pij {r ^ r + f) and p = pij (r —f r — f). 
The drift of the centre of mass therefore satisfies 



Wo 
A 

fWo 



p ■ min 1 



(49) 
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FIG. 14: Acceptance rates (vertical axis) as a function of 
maximum displacement A (horizontal axis) for two particles 
interacting with the Lennard-Jones-esque potential discussed 
in Section [nil Simulations in Section [Till were performed using 
A = O.llcr. The solid line shows the ratio of actual cluster 
displacements to the displacement predicted on the grounds of 
cluster-move proposal rate, -Ractuai/i?predictcd- The deviation 
from unity is due to the superdetailed balance factor. This 
ratio is close to unity for the displacement scale we use, in- 
dicating that little unwanted suppression of collective motion 
occurs. We show also the fraction of moves for which relative 
particle displacements are accepted (tIspm, dot-dashed line), 
and the fraction of moves for which the particles are displaced 
collectively (Aspm, curved dotted line). 



When p is small (corresponding to a small change in 
energy upon the proposed move), the collective diffusion 
of the dimer is dominated by single-particle moves, and 
hence by the first term in Equation ( 49 ) . When p ap- 



proaches 1 (corresponding to a large change in bond en- 
ergy upon the proposed move), collective diffusion in- 
duced by single-particle moves is strongly suppressed. 
Within the VMMC algorithm this suppression of dif- 
fusion is in principle compensated by explicit collective 
moves of i and j, described by the second line of Equa- 
tion ( 49 1 , provided that / and s are chosen accord- 
ingly. However, the requirement of superdetailed bal- 
ance suppresses the rate of collective moves by a fac- 
tor min {l,p/p). If this factor deviates significantly from 
unity it may by necessary to adjust the frequency / or 
displacement scale s of collective moves to restore the de- 
sired 'physical' diffusion properties of the cluster centre 
of mass. 

For the Lennard- Jones system discussed in Section III 
we may determine the unwanted suppression of cluster 
diffusion, due to the superdetailed balance factor, for two 
particles placed initially at a separation such that their 
energy of interaction is as favourable as possible. We pick 
one particle and propose a displacement of that particle. 
We draw proposed displacements f uniformly from the 



interval [0,A]. We compute the change in pair energy 
resulting from this proposed displacement, and form a 
link between particles with the virtual-move probability 
p(r —^ r+f) . If a link does not form we accept the move of 
one particle with respect to the other. If a link forms, we 
reject the proposed displacement. We instead increment 
by the 'predicted' squared displacement, ^predicted- 
this is the displacement of the dimer expected if every 
linking event results in collective motion. We then com- 
pute the reverse linking weight, p = p{r r — f), and 
with probability min(l,p/p) increment by the 'actual' 
squared displacement i^^ctuar ^'^^ * given A we pe rfor m 
this procedure 5 x 10^ times 



We show in Figure 14 



as 

a function of A the ratio of actual to predicted cluster 
displacements for moves in which a link forms between 
the two particles. We show also the fraction of moves 
in which particles move singly Aspm (no link forms) and 
collectively ^cm (a link forms and the ratio p/p permits 
movement of the dimer). Simulations in Section [lll| were 
performed using a basic translational displacement scale 
of A = O.llcr. In this regime the requirement of bal- 
ance causes little unwanted suppression of collective mo- 
tion. It should be noted that p/p is in general very close 
to unity for subsequent links made within large clusters: 
because virtual moves contains a rotational component, 
particles recruited iteratively to the pseudocluster tend 
to move larger distances (in both forward and reverse 
directions) as the iteration progresses. 

The upper limit of A is governed by the requirement 
that motion be properly 'diffusive' on length scales set by 
typical particle separations. It is likely that the utility of 
VMMC is greatest when typical particle separations are 
large compared with particle diameters, which in turn are 
large compared with potential well widths. This is the 
case for models of self-assembling proteins, whose real-life 
counterparts are typically present at low concentrations. 
It should also be noted that the effective collective mo- 
tion induced by moves of single particles produces clus- 
ter diffusivities akin to those of Brownian dynamics. If 
collective diffusion within VMMC is assigned a different 
scaling (e.g. in order to satisfy Stokes' law), it should 
be verified that the two mechanisms of cluster diffusion 
do not 'compete'. For a given model, careful considera- 
tion should be given to the choice of displacement scales 
(for both real and virtual moves) and attempted versus 
actual cluster displacements, in order to ensure that the 
desired motion of clusters is approximated. 

The generalization to many particles of the argument 
presented in this section is not straightforward. The sim- 



ple comparison performed in Section III indicates that 



for our chosen system VMMC represents a reasonable 
approximation of Brownian dynamics. The utility of 
VMMC in other cases must be assessed by testing against 
established methods, such as theoretical results (e.g. the 
solutions to Langevin equations, for sufficiently simple 
models) or Newtonian simulation. 
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